////////////////////////////TCF/////////////////////////////
// TPC Coherent Fit (TCF)                                 //
// Author: Xianguo Lu, lu.xianguo@gmail.com               //
// Reference: Xianguo Lu, "Exploring the performance      //
// limits of the ALICE Time Projection Chamber and        //
// Transition Radiation Detector for measuring identified //
// hadron production at the LHC", dissertation at the     //
// University of Heidelberg (2013).                       //
////////////////////////////////////////////////////////////

#include "dEdxGausN.h"

ClassImp(dEdxGausN);

Int_t dEdxGausN::fNtype = -999;
Double_t dEdxGausN::fMinResolution = -999;
Double_t dEdxGausN::fMaxResolution = -999;
Double_t dEdxGausN::fIniResolution = -999;

Int_t dEdxGausN::fNparType = 0;
BaseUtils::FFunc dEdxGausN::fGausFunc = 0x0;

Int_t dEdxGausN::fBBopt = -1;
Bool_t dEdxGausN::fkEnergyLossCorrection = kFALSE;
Bool_t dEdxGausN::fkSignalShift = kFALSE;
BaseUtils::FFunc dEdxGausN::fBBFunc = 0x0;
Int_t dEdxGausN::fNparBB = 0;
Double_t * dEdxGausN::fIniBBpar = 0x0;

//==================================================================================
//                               initialization
//==================================================================================
void dEdxGausN::Ini()
{
  //
  //only judged by "pwd"
  //

  const TString pwd = gSystem->pwd();

  //==================================================================================
  //                              first level from "pwd" 
  //==================================================================================
  
  //================== fBBopt
  
  fBBopt = 0;
  if(pwd.Contains("mBB"))    fBBopt += kmBB;
  if(pwd.Contains("BBXing"))    fBBopt += kBBXing;
  if(pwd.Contains("aleph"))  fBBopt += kRawALEPH;
  if(pwd.Contains("lund"))   fBBopt += kLund;

  fkEnergyLossCorrection = kFALSE;
  if(pwd.Contains("EnergyLossCorrection")){
    fkEnergyLossCorrection = kTRUE;
    printf("EnergyLossCorrection no need to be used!! %s\n", pwd.Data()); exit(1);
  }

  fkSignalShift = kFALSE;
  if(pwd.Contains("SignalShift")) fkSignalShift = kTRUE; 

  //==================================================================================
  //                              second level from internal control
  //==================================================================================

  //================== fNtype, MinResolution, MaxResolution, IniResolution

  fNtype = 4;
  if(pwd.Contains("NTYPE3"))
    fNtype = 3;

  fMinResolution = 0.03;
  fMaxResolution = 0.10;
  fIniResolution = 0.052;

  //================== fNparType, fGausFunc
  //number of parameter per type
  fNparType = 3;
  fGausFunc = BaseUtils::ResolutionGaus;

  //================== fBBFunc, fNparBB, fIniBBpar

  fNparBB = NparBasicBB();
  if(fkEnergyLossCorrection ) fNparBB++;
  if(fkSignalShift) fNparBB+= NparSignalShift();//pi,K/p
  
  fIniBBpar = new Double_t[fNparBB];
  for(Int_t ii=0; ii<fNparBB; ii++) fIniBBpar[ii]=0;

  if(fBBopt==kRawALEPH){
    fBBFunc = RawAleph;
    const Double_t tmp[]={15, 3, log(0.03), 3, 1};
    memcpy(fIniBBpar, tmp, sizeof(Double_t)*NparBasicBB());   
  }
  else if(fBBopt==kBBXing){
    /*
gSystem->Load("lib/libBaseUtils.so");
gSystem->Load("lib/libTPCHist.so");
gSystem->Load("lib/libdEdxGausN.so");
const Double_t tmp[]={40, 2,  0.2, -0.2, 80, 0,0,0,0,0,0,0,0,0};  
dEdxGausN::FindCrossingBetheBlochModified(tmp, dEdxGausN::Mkaon(), dEdxGausN::Mpion())
dEdxGausN::FindCrossingBetheBlochModified(tmp, dEdxGausN::Mkaon(), dEdxGausN::Mproton())
const Double_t tmp2[]={40, 2,  1.0492283202, 2.4938481414288, 80, 0,0,0,0,0,0,0,0,0};
TF1 *ff=new TF1("ff",dEdxGausN::BetheBlochXing, 0.1, 1e5,5)                                                                                                                              
ff->SetParameters(tmp2)
ff->SetLineColor(kBlue); ff->Draw()
ff->SetMinimum(0); ff->SetMaximum(200); gPad->SetLogx()

TF1 *f2=new TF1("f2",dEdxGausN::BetheBlochModified, 0.1, 1e5,5)
f2->SetParameters(tmp)
f2->Draw("same")
     */

    fBBFunc = BetheBlochXing;
    //const Double_t tmp[]={40, 2,  0.2, -0.2, 80, 0,0,0,0,0,0,0,0,0};
    //new p2, p3 as kaon-pion, kaon-proton crossing
    const Double_t tmp[]={40, 2,  1.0492283202, 2.4938481414288, 80};
    memcpy(fIniBBpar, tmp, sizeof(Double_t)*NparBasicBB());
  }
  else if(fBBopt==kmBB){
    fBBFunc = BetheBlochModified;
    const Double_t tmp[]={40, 2,  0.2, -0.2, 80};
    memcpy(fIniBBpar, tmp, sizeof(Double_t)*NparBasicBB());
  }
  else if(fBBopt==kLund){
    fBBFunc = Lund;
    const Double_t tmp[]={34, 0.25, 4, 1.3, 80};
    memcpy(fIniBBpar, tmp, sizeof(Double_t)*NparBasicBB());   
  }
  else{
    printf("dEdxGausN::Ini no BBopt specified wrong pwd! %s\n",pwd.Data());
    //exit(1);
  }


  Print();
}

//==================================================================================
//                                Parameters
//==================================================================================

void dEdxGausN::ResetIfixed(Bool_t pfix[])
{
  for(Int_t ipar = 0; ipar<Npar(); ipar++){
    pfix[ipar] = kFALSE;
  }
}

void dEdxGausN::SetMean(const Double_t xx, const Double_t bbpar[], Double_t pregauspar[], Int_t *neworder, const Bool_t kprint)
{
  //
  //bbpar protected by "const"
  //
  //=============== set mean ==================

  //printf("longtest setmean %e %e %e %e %e\n", bbpar[0], bbpar[1], bbpar[2], bbpar[3], bbpar[4]);

  Double_t rawys[Ntype()];
  for(Int_t itype=0; itype<Ntype(); itype++){
    rawys[itype] = GetMean(xx, Mass(itype), bbpar, kprint);
    //printf("longtest %d %e\n", itype, rawys[itype]);
  }

  Int_t id[Ntype()];
  TMath::Sort(Ntype(), rawys, id, kFALSE);

  Double_t sortedys[Ntype()];
  for(Int_t ith=0; ith<Ntype(); ith++){
    sortedys[ith] = rawys[id[ith]];
    //to repulsive parametrization of [+1]
    if(ith){
      sortedys[ith] -= rawys[id[ith-1]];
    }
  }

  for(Int_t itype=0; itype<Ntype(); itype++){
    pregauspar[itype*NparType()+1] = sortedys[itype];
  }

  if(neworder){
    //update neworder
    memcpy(neworder, id, sizeof(id));
  }
}

void dEdxGausN::ToGausPar(const Double_t *prepar, Double_t *gauspar, const Double_t *prelowerr, const Double_t *prehigherr, Double_t *gauslowerr, Double_t *gaushigherr)
{
  //dummy paramters has to be set as "fixed" in TMinuit
  //no change is made to input pre*. Protected by "const"

  //============================================================
  //prepar[NparType()*itype + 0 ]: fraction
  //============================================================

  Double_t yields[Ntype()];
  for(Int_t itype=0; itype<Ntype(); itype++){
    //keep abs, this is the internal fit par used in TMinuit.
    //i.e. even fPar[0]<0, TMinuit still see a's > 0
    yields[itype] = TMath::Abs( prepar[ itype*NparType() ] );
  }

  //============================================================
  //prepar[NparType()*itype + 1 ]: pre-mean p1, mean_i = \sum_0^i-1 abs(p1)
  //============================================================

  Double_t means[Ntype()];
  Double_t mm = 0;
  for(Int_t itype = 0; itype < Ntype(); itype++){
    //all positive
    // in this way the two peaks can not overlap and the minimum would not be the averaged peak
    mm += TMath::Abs(prepar[itype*NparType()+1]);
    means[itype] = mm;
  }

  //============================================================
  //prepar[NparType()*itype + 2 ]: always resolution
  //============================================================

  Double_t widths[Ntype()];
  for(Int_t itype=0; itype<Ntype(); itype++){
    widths[itype] = prepar[ itype*NparType()+2]; 
  }

  //============================================================
  //directory copy of other parameters
  //============================================================

  for(Int_t itype = 0; itype < Ntype(); itype++){
    gauspar[itype*NparType() + 0]=yields[itype];
    gauspar[itype*NparType() + 1]=means[itype];
    gauspar[itype*NparType() + 2]=widths[itype];
  }

  if(!prelowerr || !gauslowerr)
    return;

  SetGausErr(gauslowerr, prelowerr, prepar);
  SetGausErr(gaushigherr, prehigherr, prepar);
}

void dEdxGausN::SetGausErr(Double_t * gauserr, const Double_t * preerr, const Double_t * prepar)
{
  Double_t eM2=0;//mean error ^2
  for(Int_t itype = 0; itype < Ntype(); itype++){
    eM2 += TMath::Power(preerr[1],2);
    gauserr[itype*NparType()+1] = TMath::Sqrt(eM2);

    for(Int_t ipar = 0; ipar<NparType(); ipar++){
      if(ipar==1){
        continue;
      }

      const Int_t id = itype*NparType()+ipar;
      gauserr[id] = preerr[id];
    }
  }
}

//==================================================================================
//                             functions
//==================================================================================
Double_t dEdxGausN::ALICEBetheBlochSolid(const Double_t bg)
{
  //full credit: ALICE experiment at CERN
  //return AliExternalTrackParam::BetheBlochSolid(bg);

  const Double_t kp0=2.33;
  const Double_t kp1=0.20;
  const Double_t kp2=3.00;
  const Double_t kp3=173e-9;
  const Double_t kp4=0.49848;

  const Double_t mK  = 0.307075e-3; // [GeV*cm^2/g]                                                                                                    
  const Double_t me  = 0.511e-3;    // [GeV/c^2]                                                                                                       
  const Double_t rho = kp0;
  const Double_t x0  = kp1*2.303;
  const Double_t x1  = kp2*2.303;
  const Double_t mI  = kp3;
  const Double_t mZA = kp4;
  const Double_t bg2 = bg*bg;
  const Double_t maxT= 2*me*bg2;    // neglecting the electron mass                                                                                    

  //*** Density effect                                                                                                                                 
  Double_t d2=0.;
  const Double_t x=TMath::Log(bg);
  const Double_t lhwI=TMath::Log(28.816*1e-9*TMath::Sqrt(rho*mZA)/mI);
  if (x > x1) {
    d2 = lhwI + x - 0.5;
  } else if (x > x0) {
    const Double_t r=(x1-x)/(x1-x0);
    d2 = lhwI + x - 0.5 + (0.5 - lhwI - x0)*r*r*r;
  }

  return mK*mZA*(1+bg2)/bg2*
    (0.5*TMath::Log(2*me*bg2*maxT/(mI*mI)) - bg2/(1+bg2) - d2);

  /*
Double_t AliExternalTrackParam::BetheBlochSolid(Double_t bg) {
  //------------------------------------------------------------------                                                                                 
  // This is an approximation of the Bethe-Bloch formula,                                                                                              
  // reasonable for solid materials.                                                                                                                   
  // All the parameters are, in fact, for Si.                                                                                                          
  // The returned value is in [GeV/(g/cm^2)]                                                                                                           
  //------------------------------------------------------------------                                                                                 

  return BetheBlochGeant(bg);
}

static Double_t BetheBlochGeant(Double_t bg,
                                  Double_t kp0=2.33,
                                  Double_t kp1=0.20,
                                  Double_t kp2=3.00,
                                  Double_t kp3=173e-9,
                                  Double_t kp4=0.49848
                                  );



Double_t AliExternalTrackParam::BetheBlochGeant(Double_t bg,
         Double_t kp0,
         Double_t kp1,
         Double_t kp2,
         Double_t kp3,
         Double_t kp4) {
  //                                                                                                                                                   
  // This is the parameterization of the Bethe-Bloch formula inspired by Geant.                                                                        
  //                                                                                                                                                   
  // bg  - beta*gamma                                                                                                                                  
  // kp0 - density [g/cm^3]                                                                                                                            
  // kp1 - density effect first junction point                                                                                                         
  // kp2 - density effect second junction point                                                                                                        
  // kp3 - mean excitation energy [GeV]                                                                                                                
  // kp4 - mean Z/A                                                                                                                                    
  //                                                                                                                                                   
  // The default values for the kp* parameters are for silicon.                                                                                        
  // The returned value is in [GeV/(g/cm^2)].                                                                                                          
  //                                                                                                                                                   

  const Double_t mK  = 0.307075e-3; // [GeV*cm^2/g]                                                                                                    
  const Double_t me  = 0.511e-3;    // [GeV/c^2]                                                                                                       
  const Double_t rho = kp0;
  const Double_t x0  = kp1*2.303;
  const Double_t x1  = kp2*2.303;
  const Double_t mI  = kp3;
  const Double_t mZA = kp4;
  const Double_t bg2 = bg*bg;
  const Double_t maxT= 2*me*bg2;    // neglecting the electron mass                                                                                    

  //*** Density effect                                                                                                                                 
  Double_t d2=0.;
  const Double_t x=TMath::Log(bg);
  const Double_t lhwI=TMath::Log(28.816*1e-9*TMath::Sqrt(rho*mZA)/mI);
  if (x > x1) {
    d2 = lhwI + x - 0.5;
  } else if (x > x0) {
    const Double_t r=(x1-x)/(x1-x0);
    d2 = lhwI + x - 0.5 + (0.5 - lhwI - x0)*r*r*r;
  }

  return mK*mZA*(1+bg2)/bg2*
         (0.5*TMath::Log(2*me*bg2*maxT/(mI*mI)) - bg2/(1+bg2) - d2);
}


   */

}

Int_t dEdxGausN::NparSignalShift()
{
  return 3;
}

Double_t dEdxGausN::GetMean(const Double_t xx, const Double_t mass,  const Double_t *bbpar, const Bool_t kprint)
{
  Int_t iter = NparBasicBB();

  //----------------

  Double_t pp = TMath::Power(10,xx); //xx=log10(p)

  if(fkEnergyLossCorrection){
    const Double_t material = TMath::Abs(bbpar[iter]); iter++;
    pp = TMath::Sqrt( TMath::Power(TMath::Sqrt(pp*pp+mass*mass)-material*ALICEBetheBlochSolid(pp/mass),2) - TMath::Power(mass,2) );
  }

  //----------------

  const Double_t bg = pp/mass;

  //to speed up
  /*
  if(kprint){
    printf("getmean %e %e %e --- %e %e %e %e %e\n", TMath::Power(10,xx), pp, bg, bbpar[0], bbpar[1], bbpar[2], bbpar[3], bbpar[4]);
  }
  */

  //with +EPSILON to protect mean = 0 -> resolution -> inf
  //Max() to reject negative mean

  Double_t value = BBFunc()(&bg, bbpar);
  if(value < EPSILON)
    value = EPSILON;

  if(fkSignalShift){
    Double_t ss = 0;
    //pi
    if(     TMath::Abs(mass-Mpion())  <EPSILON){
      ss = bbpar[iter]; 
    }
    //K
    else if(TMath::Abs(mass-Mkaon())  <EPSILON){
      ss = bbpar[iter+1];
    }
    //p
    else if(TMath::Abs(mass-Mproton())<EPSILON){
      ss = bbpar[iter+2];
    }

    value *= 1 + ss;
    if(value < EPSILON)
      value = EPSILON;
  }

  if(value>TPCHist::Y1()){
    return TPCHist::Y1();
  }
  else{
    return value;
  }
}

Double_t dEdxGausN::Func(const Double_t *xx, const Double_t *inraw)
{
  Double_t par[Npar()];
  ToGausPar(inraw, par);

  return BaseUtils::GausN(xx, par, Ntype(), GausFunc(), NparType());
}

void dEdxGausN::SetBetheBlochXingKs(Double_t kk[], const Double_t a1, const Double_t pAB, const Double_t mA, const Double_t mB)
{
  const Double_t eA = TMath::Sqrt(mA*mA+pAB*pAB);
  const Double_t eB = TMath::Sqrt(mB*mB+pAB*pAB);

  const Double_t tt = BaseUtils::ProtectedPow(eB/eA, a1);

  kk[0] = 1-tt;

  //printf("test0 pAB %e mA %e mB %e eA %e eB %e a1 %e tt %e kk[0] %e\n", pAB, mA, mB, eA, eB, a1, tt, kk[0]);

  kk[1] = TMath::Log(pAB/mA) - TMath::Log(pAB/mB)*tt;
  kk[2] = TMath::Log(pAB/eA) - TMath::Log(pAB/eB)*tt;
}

void dEdxGausN::Geta2a3FromKaonCrossings(Double_t &a2, Double_t &a3, const Double_t a1, const Double_t croPi, const Double_t croProton)
{
  /*//test
gSystem->Load("lib/libBaseUtils.so");
gSystem->Load("lib/libTPCHist.so");
gSystem->Load("lib/libdEdxGausN.so");
double par[]={  3.43178e+01, 2.40223e+00,  1.90810e-01, 2.70186e-01, 7.40610e+01};
double cPi = dEdxGausN::FindCrossingBetheBlochModified(par, dEdxGausN::Mkaon(), dEdxGausN::Mpion())
double cProton = dEdxGausN::FindCrossingBetheBlochModified(par, dEdxGausN::Mkaon(), dEdxGausN::Mproton())
double a2=-999, a3=-999
dEdxGausN::Geta2a3FromKaonCrossings(a2,a3,par[1], cPi, cProton)
double as[]={par[0], par[1], a2, a3, par[4]}
double bg1=cPi/dEdxGausN::Mkaon()
double bg2=cPi/dEdxGausN::Mpion()
double y1=dEdxGausN::BetheBlochModified(&bg1, as)
double y2=dEdxGausN::BetheBlochModified(&bg2, as)
y1-y2 //should be 0
par[2]-a2 //should be 0
par[3]-a3 //should be 0

double pp[]={par[0], par[1], cPi, cProton, par[4]}
TF1 *ff=new TF1("ff",dEdxGausN::BetheBlochXing, 0.1, 1e5,5)
ff->SetParameters(pp)
ff->Draw()
   */

  Double_t K1[]={-999, -999, -999};
  Double_t K2[]={-999, -999, -999};
  SetBetheBlochXingKs(K1, a1, croPi,     Mkaon(), Mpion());
  SetBetheBlochXingKs(K2, a1, croProton, Mkaon(), Mproton());

  const Bool_t ksolve = BaseUtils::SolveEq(a2, a3, K1, K2);
  if(!ksolve){
    printf("dEdxGausN::Geta2a3FromKaonCrossings SolveEq error! K1 %e %e %e K2 %e %e %e a1 %e croPi %e croProton %e a2 %e a3 %e\n", K1[0], K1[1], K1[2], K2[0], K2[1], K2[2], a1, croPi, croProton, a2, a3); exit(1);
  }
}

Double_t dEdxGausN::BetheBlochXing(const Double_t *xx, const Double_t *rawpar)
{
  //even though this parametrization is more physics, it only introduces non-linear difficulty in minimization, otherwise no other difference
  //rawpar: r0, r1, r2,                r3,               r4
  //      = a0, a1, K-pi-crossing mom, K-p-crossing mom, a4
  //par   : a0, a1, a2,                a3,               a4
  //a1, a2, a3 from r1, r2, r3, where a1=r1

  Double_t a2 =-999, a3 = -999;
  Geta2a3FromKaonCrossings(a2, a3, rawpar[1], rawpar[2], rawpar[3]);

  const Double_t par[]={rawpar[0], rawpar[1], a2, a3, rawpar[4]};
  return BetheBlochModified(xx, par);
}

Double_t dEdxGausN::BetheBlochModified(const Double_t *xx, const Double_t *par)
{
  //===================
  //RawALEPH has numerical difficulty in ln(p2+1/pow(bg, p4))
  //use flat dedx for electron. deviation from this at high bg is not inportant because anyway there the statistics is very limited
  //npar = 5;

  const Double_t bg = xx[0];
  const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);

  const Double_t p0 = par[0];//overall prefactor
  const Double_t p1 = par[1];//beta exponent ~ 2
  const Double_t p2 = par[2];
  const Double_t p3 = par[3];
  const Double_t p4 = par[4];//plateau ~ 80

  //EPSILON added to protect 0 denominator!
  const Double_t bethebloch = p0 *( 1 + p2*TMath::Log(bg) + p3*TMath::Log(beta) )/( BaseUtils::ProtectedPow(beta, p1)+EPSILON) ;

  //tp = 100, electron momentum = 51 MeV, pion momentum = 14 GeV
  const Double_t tp = 10;

  Double_t value = 0;
  if(bg<tp){
    value = bethebloch;
  }
  else{
    value = TMath::Min(p4, bethebloch);
  }

  return value;
}

//const Double_t x0 = dEdxGausN::FindCrossing(fPars, , dEdxGausN::kPION,   dEdxGausN::kKAON);

Double_t dEdxGausN::FindCrossingLog10P(const Double_t* par, BaseUtils::FFunc func, const Double_t m1, const Double_t m2)
{
  //return log10(beta gamma)

  TH1D *hh=new TH1D("tmpFindCrossing","",5000, TMath::Log10(0.5), TMath::Log10(5));

  for(Int_t ii=1; ii<=hh->GetNbinsX(); ii++){
    const Double_t xx = hh->GetBinCenter(ii);
    const Double_t bg1 = TMath::Power(10,xx)/m1;
    const Double_t bg2 = TMath::Power(10,xx)/m2;

    const Double_t y1 = func(&bg1, par);
    const Double_t y2 = func(&bg2, par);
    hh->SetBinContent(ii, -TMath::Abs(y1-y2));
  }

  const Double_t mx = hh->GetBinCenter(hh->GetMaximumBin());

  delete hh;

  return mx;
}

Double_t dEdxGausN::FindCrossingBetheBlochModified(const Double_t *par, const Double_t m1, const Double_t m2)
{
  const Double_t mx = FindCrossingLog10P(par, BetheBlochModified, m1, m2);
  return TMath::Power(10,mx);
}

Double_t dEdxGausN::RawAleph(const Double_t *xx, const Double_t *par)
{
  const Double_t bg = xx[0];
  const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);

  const Double_t p0 = TMath::Abs(par[0]);
  const Double_t p1 = TMath::Abs(par[1]);
  const Double_t p2 = par[2];//TMath::Abs(par[2]);//still first have to be positive
  const Double_t p3 = TMath::Abs(par[3]);
  const Double_t p4 = TMath::Abs(par[4]);//has to be positive so that in bg-> inf there is no bg dependence anymore

  //return p0*((p1-TMath::Log(p2+1/BaseUtils::ProtectedPow(bg,p4)))/BaseUtils::ProtectedPow(beta,p3) - 1 );
  //return p0*( (p1- TMath::Log(p2+EPSILON)- TMath::Log(1+1/(BaseUtils::ProtectedPow(bg,p4)*p2+EPSILON))  )/( BaseUtils::ProtectedPow(beta,p3)+EPSILON ) - 1 );
  //log(p2) -> p2 because p2 sometimes ~EPSILON

  //truncate overflow so that the fraction is suppressed as 0 for non-existing peak
  //should well above non-zero signal

  const Double_t value = p0*( (p1- p2 - TMath::Log(1+1/(BaseUtils::ProtectedPow(bg,p4)*BaseUtils::ProtectedExp(p2)+EPSILON))  )/( BaseUtils::ProtectedPow(beta,p3)+EPSILON ) - 1 );

  return value;

  //p2 can be very small ~ 1e-9, do not separate p2 out! -> absorb p2 into p1, remember to change initial value of p1
  //return p0*( (p1- TMath::Log(1+1/(BaseUtils::ProtectedPow(bg,p4)*p2+EPSILON))  )/( BaseUtils::ProtectedPow(beta,p3)+EPSILON ) - 1 );
  //change p0 as plateau
  //return p0/(p1-1)*( (p1- TMath::Log(1+1/(BaseUtils::ProtectedPow(bg,p4)*p2+EPSILON))  )/( BaseUtils::ProtectedPow(beta,p3)+EPSILON ) - 1 );
  //redefine p1 so that p1>1 --> p1>0
  //return p0/p1*( (p1 + 1 - TMath::Log(1+1/(BaseUtils::ProtectedPow(bg,p4)*p2+EPSILON))  )/( BaseUtils::ProtectedPow(beta,p3)+EPSILON ) - 1 );
}

Double_t dEdxGausN::Lund(const Double_t * xx,  const Double_t * rawpar)
{
  Double_t par[100];
  memcpy(&(par[1]), rawpar, sizeof(Double_t)*5);

  const Double_t bg = xx[0];

  // bg is beta*gamma
  const Double_t beta2 = bg*bg / (1.0+bg*bg);
  const Double_t a = par[1];
  const Double_t b = par[2]*par[1];//in this way value = rawpar[0]*()
  const Double_t c = par[3];
  const Double_t e = par[4];
  const Double_t f = par[5];
  const Double_t d = BaseUtils::ProtectedExp(c*(a-f)/b);

  const Double_t powbg = BaseUtils::ProtectedPow(1.0+bg, c);

  const Double_t value = a/BaseUtils::ProtectedPow(beta2,e) + b/c*TMath::Log(powbg/(1.0 + d*powbg));

  return value;
}

//==================================================================================
//                             algorithms
//==================================================================================
void dEdxGausN::GetInvMap(const Int_t order[], Int_t invmap[])
{
  GetInvMap(Ntype(), order, invmap);
}

void dEdxGausN::GetInvMap(const Int_t nn, const Int_t order[], Int_t invmap[])
{
 for(Int_t itype=0; itype<nn; itype++){
    invmap[itype]=-999;
    for(Int_t ss=0; ss<nn; ss++){
      if(order[ss]==itype){
        invmap[itype]=ss;
        break;
      }
    }
    if( invmap[itype]<0 ){
      printf("dEdxGausN::GetInvMap fail to find invmap %d!\n", itype);exit(1);
    }
  }
}
void dEdxGausN::Swap(const Double_t xx, const Int_t targetip, Double_t pars[], const Int_t oldid[], const Int_t newid[])
{
  Double_t oldpar[Npar()];
  memcpy(oldpar, pars, sizeof(oldpar));
  
  for(Int_t newith=0; newith<Ntype(); newith++){
    Int_t oldith = -999;
    for(Int_t ss=0; ss<Ntype(); ss++){
      if(oldid[ss]==newid[newith]){  //the newith-th par now = the ss-th par before
        oldith=ss;
        break;
      }
    }
    if(oldith<0){
      printf("dEdxGausN::Swap old id not found!! error!!\n");exit(1);
    }
    
    if(newith!=oldith){
      printf("dEdxGausN::Swap reordering %d at xx %e: type %d old %d new %d\n", targetip, xx, newid[newith], oldith, newith);
     
      pars[NparType()* newith + targetip] = oldpar[NparType()* oldith +targetip];
    }
  }
  
}

//==================================================================================
//                               kinematics
//==================================================================================
Double_t dEdxGausN::Mass(const Int_t id)
{
  switch(id){
  case kPROTON:
    return Mproton();
  case kPION:
    return Mpion();
  case kKAON:
    return Mkaon();
  case kELECTRON:
    return Melectron();
  default:
    printf("dEdxGausN::Mass wrong id %d\n", id); exit(1);
  }

  return -999;
}

const Char_t * dEdxGausN::ParticleName(const Int_t id)
{
 switch(id){
  case kPROTON:
    return "Proton";
  case kPION:
    return "Pion";
  case kKAON:
    return "Kaon";
  case kELECTRON:
    return "Electron";
  default:
    printf("dEdxGausN::ParticleName wrong id %d\n", id); exit(1);
  }

 return "xxx";
}

//==================================================================================
//                            histogram
//==================================================================================

//==================================================================================
//                            outputs
//==================================================================================

void dEdxGausN::Print()
{
  printf("\ndEdxGausN::Print()\n");

  printf("important controls:\n");

  printf("fBBopt %d\n", fBBopt);
  printf("fkEnergyLossCorrection %d\n", fkEnergyLossCorrection);
  printf("fkSignalShift %d\n", fkSignalShift);
  printf("NparSignalShift() %d\n", NparSignalShift());
  printf("\n");

  printf("Ntype(): %d\n", Ntype());
  printf("MinResolution(): %e\n", MinResolution());
  printf("MaxResolution(): %e\n", MaxResolution());
  if(MinResolution()<0 || MaxResolution()<0){
    printf("dEdxGausN::Print min/max res < 0 !\n"); exit(1);
  }

  printf("IniResolution(): %e\n", IniResolution());

  printf("NparType(): %d\n", NparType());

  printf("NparBB(): %d\n", NparBB());
  printf("TPCHist::Y1(): %e\n", TPCHist::Y1());
  printf("\n");
}

Int_t dEdxGausN::SetIfixAllBB(Bool_t pfix[], const Bool_t kfix)
{
  const Int_t nfix = NparBB();
  Int_t nflip = 0;
  for(Int_t ii=0; ii< nfix; ii++){
    if(pfix[ii]!=kfix){
      pfix[ii] = kfix;
      nflip++;
    }
  }
  if(kfix)
    return nflip;
  else
    return nflip*(-1);
}

Int_t dEdxGausN::SetIfixExtraBB(Bool_t pfix[], const Bool_t kfix)
{
  const Int_t nfix = NparBB()-NparBasicBB();
  Int_t nflip = 0;
  for(Int_t ii=0; ii< nfix; ii++){
    const Int_t dd = NparBB()-1-ii;
    if(pfix[dd]!=kfix){
      pfix[dd] = kfix;
      nflip++;
    }
  }

  if(kfix)
    return nflip;
  else
    return nflip*(-1);
}

TH1D *dEdxGausN::dEdxFit(TGraphErrors * gr, const TString hname, const Bool_t kfailexit, const Double_t *inipar, const Double_t *lows, const Double_t *highs, const Double_t maxchi2)
{
  printf("dEdxGausN::dEdxFit %s start ...\n", hname.Data());

  Double_t par[100];
  if(!inipar){
    inipar = fIniBBpar;
  }
  memcpy(par, inipar, sizeof(Double_t)* NparBB());
  Double_t err[100], chi[100];

  Bool_t pfix[NparBB()];
  SetIfixAllBB(pfix, kFALSE);

  if(fkEnergyLossCorrection || fkSignalShift){
    SetIfixExtraBB(pfix, kTRUE);
  }

  const Int_t kfail = BaseUtils::ChisquareFit(gr->GetN(), gr->GetX(), gr->GetY(), 0x0, gr->GetEY(), BBFunc(), NparBB(), par,err,chi, pfix, 0x0, lows, highs); 
 
  if(kfail || chi[0]/chi[1] > maxchi2){
    printf("dEdxGausN::dEdxFit kfail %d %s or bad chi2 %e %e, %e %e\n\n", kfail,kfailexit?"exit!!!":"", chi[0], chi[1], chi[0]/chi[1], maxchi2);
    if(kfailexit){
      printf("exit!!!\n");
      exit(1);
    }
  }

  /*
  if(fBBopt==kRawALEPH){
    if(par[0]>100){
      printf("dEdxGausN::dEdxFit special check for p0 aleph, larger than 100, fall back to inipar! %e -- %e\n", par[0], inipar[0]);
      memcpy(par, inipar, sizeof(Double_t)* NparBB());
    }
  }
  */

  TH1D *hfit = BaseUtils::GetHfit(hname, BBFunc(), par, 0.1, 4e4, kTRUE);
  TString htit;
  for(Int_t ii=0; ii<NparBB(); ii++){
    htit+=Form("%e , ",par[ii]);
  }
  hfit->SetTitle(htit);
  return hfit;
}

void dEdxGausN::DrawTF1(const Int_t id, const TString tag,  const Double_t pars[], const Int_t col, TList *ll, FILE *log)
{
  Char_t buff[2000];
  sprintf(buff, "dEdxGasuN::DrawTF1 longtest %s %d %e %e %e\n", tag.Data(), id, pars[0], pars[1], pars[2]);
  if(log){
    fprintf(log, buff,"");
  }
  else{
    printf(buff,"");
  }

  TH1D * ff = BaseUtils::GetHfit(Form("%sf%d", tag.Data(), id), GausFunc(), pars, 0, 200);
  ff->SetLineColor(col);
  //ff->SetDrawOption("same hist C");
  ll->Add(ff);
}

void dEdxGausN::SavePar(TTreeSRedirector * stream, const TString tn,  Double_t *totalyield,  Double_t *yields, Double_t *pars,  Double_t *lowerrs,  Double_t *higherrs, const Int_t order[], Double_t *outpar, Double_t *outlowerr, Double_t *outhigherr, Int_t *orderinv)
{
  //
  //totalyield, pars, lowerrs, higherrs should not be changed
  //due to streamer, no protection is made using "const"
  //

  (*stream)<<tn<<
    "totn="<<*totalyield;

  //========================================================================
  //========================================================================

  for(Int_t itype=0; itype<Ntype(); itype++){
    for(Int_t ipar = 0; ipar<NparType(); ipar++){
      (*stream)<<tn<<
        Form("rawp%d%d=",itype, ipar)<<pars[itype*NparType()+ipar]<<
        Form("rawle%d%d=",itype, ipar)<<lowerrs[itype*NparType()+ipar]<<
        Form("rawhe%d%d=",itype, ipar)<<higherrs[itype*NparType()+ipar];
    }
  }

  //========================================================================
  //                       the following only change outpar
  //========================================================================

  //outpar and outerr have to be define outside so to have a fixed memory address because the streamer does not end here. It uses the memory address after this function.
  //pars needs to be used for next fit, can not be touched!

  ToGausPar(pars, outpar, lowerrs, higherrs, outlowerr, outhigherr);

  //========================================================
  //                         save
  //========================================================

  //the inverse map of order: invmap[itype] = ith; order[ith]=itype
  Int_t invmap[Ntype()];
  GetInvMap(order, invmap);

  for(Int_t itype=0; itype<Ntype(); itype++){
    Double_t &iyield = yields[invmap[itype]];
    iyield = *totalyield*outpar[invmap[itype]*NparType()];

    (*stream)<<tn<<
      Form("y%d=",itype)<<iyield;

    for(Int_t ipar = 0; ipar<NparType(); ipar++){
      //-> the oder of streaming have to stay always the same
      //-> save in order of particle type

      //only address of outpar and outerr can be used for streamer!
      Double_t &var = outpar[invmap[itype]*NparType()+ipar];
      Double_t &lowerr = outlowerr[invmap[itype]*NparType()+ipar];
      Double_t &higherr= outhigherr[invmap[itype]*NparType()+ipar];
      
      (*stream)<<tn<<
        Form("p%d%d=",itype, ipar)<<var<<
        Form("le%d%d=",itype, ipar)<<lowerr<<
        Form("he%d%d=",itype, ipar)<<higherr;
    }
  }

  if(orderinv){
    memcpy(orderinv, invmap, sizeof(invmap));
  }
}
